Class 2

Zanim zaczniemy

0. Przydatne zasoby pomocowe w pracy z `R` - część druga

A. `tidyverse` oraz `dplyr` - “porządne dane”

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

Ekosystem tidyverse to niejako “państwo w państwie” - rodzaj dialektu R w mocny sposób ingerujący w sposób, w jaki obsługiwane są pewne stałe elementy języka. Podstawowa rola tidyverse to stworzenie wspólnego i spójnego języka, w którym można opisywać standardowe struktury danych (tibble, dplyr), automatyzować funkcje i działania (purrr) oraz standardowo wizualizować dane w oparciu o gramatykę grafiki (ggplot2).

1. Filtrowanie danych: `filter()`

Filtrowanie odbywa się za pomocą par zmienna OPERATOR warunek - podobnie jak przy filtrowaniu za pomocą indeksów. dplyr upraszcza wszystko pozwalając nam bezpośrednio odwoływać się do nazw zmiennych.

Przefiltrujmy dane tak, by wyświetlić tylko potomstwo matki nr. R187528:

dane <- read.table(file = "BTdata.csv",
                   sep = ",", header = T,
                   stringsAsFactors = T)

filter(dane, dam == "R187528")

##   tarsus     back     dam fosternest hatchdate  sex weight habitat bill_length
## 1  14.46 549.6392 R187528      A2602        45  Fem    9.7  forest      11.565
## 2  14.61 549.6629 R187528      A2602        45  Fem   10.0  forest      11.737
## 3  14.16 549.3603 R187528      A2602        45  Fem    9.4    park      11.378
## 4  14.08 552.7623 R187528      A2602        45  Fem    9.4  forest      11.261
## 5  14.08 552.5236 R187528      A2602        45  Fem    9.5    park      11.296
## 6  14.77 550.0366 R187528      A1302        45 Male   10.1    park      11.825
## 7  14.01 551.5308 R187528      A1302        45  Fem    9.7    park      11.225
## 8  15.07 549.1769 R187528      A1302        45 Male   10.1    park      12.100
##   bill_depth
## 1      0.631
## 2      0.668
## 3      0.583
## 4      0.603
## 5      0.593
## 6      0.652
## 7      0.631
## 8      0.684

Czy wszystkie młode wychowywały się w jednym gnieździe (fosternest)?

unique(filter(dane, dam == "R187528")$fosternest)

## [1] A2602 A1302
## 104 Levels: A1002 A102 A1202 A1302 A1602 A1802 A18B02 A2202 A22B02 ... H702

Wybierzmy spośród tych piskląt tylko takie, których skok przekracza 14.50 mm, lub masa przekracza, lub jest równa, 9.5 g:

filter(dane, dam == "R187528" & (tarsus > 14.50 | weight >= 9.5))

##   tarsus     back     dam fosternest hatchdate  sex weight habitat bill_length
## 1  14.46 549.6392 R187528      A2602        45  Fem    9.7  forest      11.565
## 2  14.61 549.6629 R187528      A2602        45  Fem   10.0  forest      11.737
## 3  14.08 552.5236 R187528      A2602        45  Fem    9.5    park      11.296
## 4  14.77 550.0366 R187528      A1302        45 Male   10.1    park      11.825
## 5  14.01 551.5308 R187528      A1302        45  Fem    9.7    park      11.225
## 6  15.07 549.1769 R187528      A1302        45 Male   10.1    park      12.100
##   bill_depth
## 1      0.631
## 2      0.668
## 3      0.593
## 4      0.652
## 5      0.631
## 6      0.684

Zadanie Wybierz wszystkie obserwacje, które są samicami i mają kolor płaszcza (back) o dominancie barwnej 550 nm (w zaokrągleniu w dół do liczb całkowitych; użyj floor()), oraz wykluły się (hatchdate) przed dniem 45 (15 maja). Ile jest takich obserwacji?

Spodziewany wynik

filter(dane, sex == "Fem" & floor(back) == 550 & hatchdate < 45)

##   tarsus     back     dam fosternest hatchdate sex weight habitat bill_length
## 1  14.31 550.0810 R187517       C602        44 Fem    9.5    park      11.434
## 2  15.22 550.1262 R187517       C602        44 Fem   10.1    park      12.206
## 3  14.61 550.2081 R187517       C602        44 Fem    9.5  forest      11.698
## 4  14.54 550.1580 R187517       C602        44 Fem    9.8    park      11.620
## 5  14.99 550.6216 R187517      C2202        44 Fem   10.3    park      12.042
## 6  14.69 550.4395 R187398       C602        44 Fem   10.1  forest      11.801
## 7  14.99 550.3093 R187517      C2202        44 Fem   10.4    park      12.021
##   bill_depth
## 1      0.629
## 2      0.696
## 3      0.644
## 4      0.573
## 5      0.679
## 6      0.638
## 7      0.677

nrow(filter(dane, sex == "Fem" & floor(back) == 550 & hatchdate < 45))

## [1] 7

Zadanie Zmodyfikuj powyższy filtr, by wyciągnąć samice o danej dominancie barwnej LUB te wyklute przed 15 maja. Ile jest takich obserwacji?

Spodziewany wynik

nrow(filter(dane, sex == "Fem" & (floor(back) == 550 | hatchdate < 45)))

## [1] 90

$\blacksquare$

2. Sortowanie za pomocą `arrange()`

Posortujmy dane mające dominantę barwną 550 nm (jak wyżej) oraz wyklute przed 15 maja wg długości skoku.

dane2 <- filter(dane, floor(back) == 550 & hatchdate < 45)
arrange(dane2, tarsus)

##    tarsus     back     dam fosternest hatchdate  sex weight habitat bill_length
## 1   14.16 550.4797 R187398       C602        44 Male    9.6  forest      11.295
## 2   14.27 550.5134 R046108       C402        44  UNK    9.6  forest      11.403
## 3   14.31 550.0810 R187517       C602        44  Fem    9.5    park      11.434
## 4   14.54 550.1580 R187517       C602        44  Fem    9.8    park      11.620
## 5   14.54 550.2464 R187517      C2202        44 Male    9.8  forest      11.644
## 6   14.61 550.2081 R187517       C602        44  Fem    9.5  forest      11.698
## 7   14.61 550.2466 R187517      C2202        44 Male    9.8    park      11.623
## 8   14.69 550.4395 R187398       C602        44  Fem   10.1  forest      11.801
## 9   14.69 550.3631 R187398       C602        44 Male    9.8  forest      11.719
## 10  14.99 550.6216 R187517      C2202        44  Fem   10.3    park      12.042
## 11  14.99 550.3093 R187517      C2202        44  Fem   10.4    park      12.021
## 12  15.22 550.1262 R187517       C602        44  Fem   10.1    park      12.206
##    bill_depth
## 1       0.644
## 2       0.598
## 3       0.629
## 4       0.573
## 5       0.654
## 6       0.644
## 7       0.622
## 8       0.638
## 9       0.679
## 10      0.679
## 11      0.677
## 12      0.696

Zmodyfikujmy sortowanie tak, by uporządkowanie według długości skoku odbywało się w obrębie każdej matki genetycznej (dam).

arrange(dane2, dam, tarsus)

##    tarsus     back     dam fosternest hatchdate  sex weight habitat bill_length
## 1   14.27 550.5134 R046108       C402        44  UNK    9.6  forest      11.403
## 2   14.16 550.4797 R187398       C602        44 Male    9.6  forest      11.295
## 3   14.69 550.4395 R187398       C602        44  Fem   10.1  forest      11.801
## 4   14.69 550.3631 R187398       C602        44 Male    9.8  forest      11.719
## 5   14.31 550.0810 R187517       C602        44  Fem    9.5    park      11.434
## 6   14.54 550.1580 R187517       C602        44  Fem    9.8    park      11.620
## 7   14.54 550.2464 R187517      C2202        44 Male    9.8  forest      11.644
## 8   14.61 550.2081 R187517       C602        44  Fem    9.5  forest      11.698
## 9   14.61 550.2466 R187517      C2202        44 Male    9.8    park      11.623
## 10  14.99 550.6216 R187517      C2202        44  Fem   10.3    park      12.042
## 11  14.99 550.3093 R187517      C2202        44  Fem   10.4    park      12.021
## 12  15.22 550.1262 R187517       C602        44  Fem   10.1    park      12.206
##    bill_depth
## 1       0.598
## 2       0.644
## 3       0.638
## 4       0.679
## 5       0.629
## 6       0.573
## 7       0.654
## 8       0.644
## 9       0.622
## 10      0.679
## 11      0.677
## 12      0.696

Zadanie Posortuj cały zestaw danych według wagi ciała, a następnie według gniazda fosternest. Alfabetycznie - jakie gniazda wypadają w najniższej kategorii wagowej? (Spróbuj wyświetlić tylko pierwsze wiersze posortowanego zestawu danych, reszty nie potrzebujesz - np. 10-20 pierwszych).

Spodziewany efekt

arrange(dane, weight, fosternest)[1:20,]

##    tarsus     back     dam fosternest hatchdate  sex weight habitat bill_length
## 1   13.10 549.4057 R186912      B1102        54  Fem    8.6    park      10.467
## 2   12.57 551.1584 R186908      B1602        53  Fem    8.6    park      10.071
## 3   13.10 549.2661 R187539       A102        47  Fem    8.8  forest      10.416
## 4   13.25 548.5838 R186912      B1102        54  Fem    8.8    park      10.559
## 5   12.95 550.3171 R187523      E2002        46  Fem    8.8    park      10.332
## 6   13.10 551.0305 R046109       C402        44  Fem    8.9    park      10.448
## 7   13.25 550.9402 R187557      F1902        47 Male    8.9  forest      10.600
## 8   13.55 551.3757 R187557      F2102        47  Fem    8.9    park      10.885
## 9   13.40 549.0100 R187588      F2402        47  Fem    8.9    park      10.674
## 10  13.48 548.0596 K983388      H1302        48  Fem    8.9  forest      10.783
## 11  13.25 549.0567   Fem20      H1302        48  Fem    8.9  forest      10.558
## 12  13.40 549.0138 R187527      A1802        45  Fem    9.0    park      10.772
## 13  13.02 554.1662 R187562      A2302        47  Fem    9.0  forest      10.441
## 14  13.25 550.1252 R187562       A602        47  Fem    9.0    park      10.616
## 15  13.55 549.6995 R187902      B1102        54  UNK    9.0  forest      10.850
## 16  13.86 549.4298 R187595      C2402        48  Fem    9.0  forest      11.034
## 17  13.33 548.9120 R186907      E1002        53  Fem    9.0    park      10.648
## 18  13.86 548.7900 R187531      F1702        46  Fem    9.0  forest      11.106
## 19  13.48 550.5093 R187916      G1602        49  Fem    9.0    park      10.798
## 20  13.78 549.4879 R187931      G2202        49  Fem    9.0    park      11.008
##    bill_depth
## 1       0.584
## 2       0.574
## 3       0.573
## 4       0.526
## 5       0.517
## 6       0.587
## 7       0.580
## 8       0.597
## 9       0.577
## 10      0.620
## 11      0.548
## 12      0.593
## 13      0.632
## 14      0.601
## 15      0.577
## 16      0.570
## 17      0.600
## 18      0.599
## 19      0.604
## 20      0.647

Zadanie Zmodyfikuj powyższe wyszukiwanie zmieniając cechę na skok, oraz sortując malejąco (użyj modyfikatora desc na zmiennej tarsus).

Spodziewany wynik

arrange(dane, desc(tarsus), fosternest)[1:20,]

##    tarsus     back     dam fosternest hatchdate  sex weight habitat bill_length
## 1   16.81 548.3708 R187545      D1002        47 Male   11.4    park      13.399
## 2   15.98 549.7797 R187537      A2702        47 Male   10.8    park      12.828
## 3   15.67 550.5012 R187914       G602        49 Male   10.5    park      12.535
## 4   15.64 550.5361 R187292      E1802        47 Male   10.2  forest      12.471
## 5   15.60 550.4936 R187155      F2402        47 Male   11.2    park      12.477
## 6   15.52 551.3060 R186908      B1602        53  UNK   10.4  forest      12.391
## 7   15.52 549.4172 R186908      B1702        53 Male   10.4    park      12.418
## 8   15.52 550.8471 R187030      C1302        46 Male   10.5  forest      12.403
## 9   15.52 551.7863 R187590       D202        49 Male   10.5  forest      12.394
## 10  15.45 548.3911 R187030      C1302        46  Fem   10.2  forest      12.348
## 11  15.45 549.9689 R187030      C1302        46 Male   10.6  forest      12.324
## 12  15.45 550.7465 R187927       D402        50 Male    9.8    park      12.298
## 13  15.45 551.0630 R187577      E1802        47 Male   10.4    park      12.358
## 14  15.45 550.6196 R187521      E2002        46 Male   10.4    park      12.345
## 15  15.45 549.8546 R187399       E402        47 Male   10.4    park      12.350
## 16  15.45 552.1680 R187531      F1702        46 Male   10.4  forest      12.378
## 17  15.45 549.6147 R187030      F1702        46 Male   10.6    park      12.424
## 18  15.45 548.6120 R187030      F1702        46 Male   10.3  forest      12.384
## 19  15.45 547.8379   Fem20      H1302        48 Male   10.8    park      12.389
## 20  15.37 549.5868 R187515      A1802        45 Male   10.5    park      12.289
##    bill_depth
## 1       0.797
## 2       0.690
## 3       0.711
## 4       0.709
## 5       0.696
## 6       0.693
## 7       0.674
## 8       0.697
## 9       0.725
## 10      0.675
## 11      0.706
## 12      0.671
## 13      0.682
## 14      0.624
## 15      0.669
## 16      0.669
## 17      0.675
## 18      0.696
## 19      0.728
## 20      0.633

$\blacksquare$

3. Selekcja zmiennych dzięki `select()`

select() pozwala wybierać pojedyncze kolumny z danych, albo za pomocą ich pełnej nazwy albo za pomocą modyfikatorów dopasowujących kolumny do zadanych warunków:

starts_with() - dopasowuje kolumna na podstawie początkowych znaków
ends_with() - podobnie, jak wyżej, ale dla znaków końcowych
contains() - dopasowuje na podstawie ciągu znaków zawartego gdzieś w nazwie kolumny
matches() - dopasowuje kolumny na podstawie wyrażenia regularnego (patrz linki na górze)

Wybierz z zestawu danych kolumny od 1 do 4 oraz 8 (we wszystkich przykładach wyświetlamy pierwszych 20 wierszy za pomocą head()):

head(select(dane, 1:4, 8),
     20)

##    tarsus     back     dam fosternest habitat
## 1   13.55 551.3757 R187557      F2102    park
## 2   15.07 549.0884 R187559      F1902  forest
## 3   14.99 550.1739 R187568       A602  forest
## 4   14.69 550.3067 R187518      A1302  forest
## 5   14.46 549.6392 R187528      A2602  forest
## 6   13.93 551.8693 R187945      C2302  forest
## 7   13.93 549.4878    Fem3      C1902    park
## 8   15.45 548.3911 R187030      C1302  forest
## 9   14.31 550.0810 R187517       C602    park
## 10  14.46 549.8523 R187523      B2202  forest
## 11  14.99 547.7403 R186902      B1402  forest
## 12  14.08 548.9847 R187400      B1002    park
## 13  14.77 551.1920 R187932       B502    park
## 14  13.63 549.2743 R187582      D1202  forest
## 15  16.81 548.3708 R187545      D1002    park
## 16  14.24 551.1921 R187546       D902    park
## 17  14.92 551.5066 R187590       D202  forest
## 18  14.46 549.7494 R187548       E902    park
## 19  15.37 552.3874 R187594       E302  forest
## 20  13.40 549.0100 R187588      F2402    park

Wybierz kolumny zawierające pomiary strukturalne (skok oraz dwa pomiary dzioba):

head(select(dane, tarsus, starts_with("bill")),
     20)

##    tarsus bill_length bill_depth
## 1   13.55      10.885      0.597
## 2   15.07      12.056      0.682
## 3   14.99      12.000      0.673
## 4   14.69      11.733      0.639
## 5   14.46      11.565      0.631
## 6   13.93      11.165      0.626
## 7   13.93      11.154      0.591
## 8   15.45      12.348      0.675
## 9   14.31      11.434      0.629
## 10  14.46      11.564      0.609
## 11  14.99      11.977      0.636
## 12  14.08      11.251      0.685
## 13  14.77      11.828      0.638
## 14  13.63      10.904      0.612
## 15  16.81      13.399      0.797
## 16  14.24      11.409      0.679
## 17  14.92      11.929      0.700
## 18  14.46      11.512      0.669
## 19  15.37      12.259      0.641
## 20  13.40      10.674      0.577

Przenieś kolumnę hatchdate na sam koniec tabeli:

head(select(dane, -hatchdate, everything(), hatchdate),
     20)

##    tarsus     back     dam fosternest  sex weight habitat bill_length
## 1   13.55 551.3757 R187557      F2102  Fem    8.9    park      10.885
## 2   15.07 549.0884 R187559      F1902 Male   10.5  forest      12.056
## 3   14.99 550.1739 R187568       A602 Male    9.9  forest      12.000
## 4   14.69 550.3067 R187518      A1302 Male    9.7  forest      11.733
## 5   14.46 549.6392 R187528      A2602  Fem    9.7  forest      11.565
## 6   13.93 551.8693 R187945      C2302  Fem    9.3  forest      11.165
## 7   13.93 549.4878    Fem3      C1902 Male    9.3    park      11.154
## 8   15.45 548.3911 R187030      C1302  Fem   10.2  forest      12.348
## 9   14.31 550.0810 R187517       C602  Fem    9.5    park      11.434
## 10  14.46 549.8523 R187523      B2202  Fem    9.6  forest      11.564
## 11  14.99 547.7403 R186902      B1402 Male   10.0  forest      11.977
## 12  14.08 548.9847 R187400      B1002  Fem    9.7    park      11.251
## 13  14.77 551.1920 R187932       B502 Male    9.9    park      11.828
## 14  13.63 549.2743 R187582      D1202 Male    9.2  forest      10.904
## 15  16.81 548.3708 R187545      D1002 Male   11.4    park      13.399
## 16  14.24 551.1921 R187546       D902  Fem    9.8    park      11.409
## 17  14.92 551.5066 R187590       D202  Fem   10.1  forest      11.929
## 18  14.46 549.7494 R187548       E902  Fem    9.6    park      11.512
## 19  15.37 552.3874 R187594       E302 Male   10.4  forest      12.259
## 20  13.40 549.0100 R187588      F2402  Fem    8.9    park      10.674
##    bill_depth hatchdate
## 1       0.597        47
## 2       0.682        47
## 3       0.673        47
## 4       0.639        45
## 5       0.631        45
## 6       0.626        49
## 7       0.591        47
## 8       0.675        46
## 9       0.629        44
## 10      0.609        46
## 11      0.636        52
## 12      0.685        47
## 13      0.638        49
## 14      0.612        48
## 15      0.797        47
## 16      0.679        47
## 17      0.700        49
## 18      0.669        48
## 19      0.641        49
## 20      0.577        47

Zadanie Usuń z tabeli danych kolumny z pomiarami morfometrycznymi (skok, waga, dziób).

Spodziewany wynik

head(select(dane, -tarsus, weight, -starts_with("bill")),
     20)

##        back     dam fosternest hatchdate  sex weight habitat
## 1  551.3757 R187557      F2102        47  Fem    8.9    park
## 2  549.0884 R187559      F1902        47 Male   10.5  forest
## 3  550.1739 R187568       A602        47 Male    9.9  forest
## 4  550.3067 R187518      A1302        45 Male    9.7  forest
## 5  549.6392 R187528      A2602        45  Fem    9.7  forest
## 6  551.8693 R187945      C2302        49  Fem    9.3  forest
## 7  549.4878    Fem3      C1902        47 Male    9.3    park
## 8  548.3911 R187030      C1302        46  Fem   10.2  forest
## 9  550.0810 R187517       C602        44  Fem    9.5    park
## 10 549.8523 R187523      B2202        46  Fem    9.6  forest
## 11 547.7403 R186902      B1402        52 Male   10.0  forest
## 12 548.9847 R187400      B1002        47  Fem    9.7    park
## 13 551.1920 R187932       B502        49 Male    9.9    park
## 14 549.2743 R187582      D1202        48 Male    9.2  forest
## 15 548.3708 R187545      D1002        47 Male   11.4    park
## 16 551.1921 R187546       D902        47  Fem    9.8    park
## 17 551.5066 R187590       D202        49  Fem   10.1  forest
## 18 549.7494 R187548       E902        48  Fem    9.6    park
## 19 552.3874 R187594       E302        49 Male   10.4  forest
## 20 549.0100 R187588      F2402        47  Fem    8.9    park

Zadanie Pozostaw w tabeli tylko kolumny zawierające w nazwie przynajmniej jedno “e”.

Spodziewany wynik

head(select(dane, contains("e")),
     20)

##    fosternest hatchdate  sex weight bill_length bill_depth
## 1       F2102        47  Fem    8.9      10.885      0.597
## 2       F1902        47 Male   10.5      12.056      0.682
## 3        A602        47 Male    9.9      12.000      0.673
## 4       A1302        45 Male    9.7      11.733      0.639
## 5       A2602        45  Fem    9.7      11.565      0.631
## 6       C2302        49  Fem    9.3      11.165      0.626
## 7       C1902        47 Male    9.3      11.154      0.591
## 8       C1302        46  Fem   10.2      12.348      0.675
## 9        C602        44  Fem    9.5      11.434      0.629
## 10      B2202        46  Fem    9.6      11.564      0.609
## 11      B1402        52 Male   10.0      11.977      0.636
## 12      B1002        47  Fem    9.7      11.251      0.685
## 13       B502        49 Male    9.9      11.828      0.638
## 14      D1202        48 Male    9.2      10.904      0.612
## 15      D1002        47 Male   11.4      13.399      0.797
## 16       D902        47  Fem    9.8      11.409      0.679
## 17       D202        49  Fem   10.1      11.929      0.700
## 18       E902        48  Fem    9.6      11.512      0.669
## 19       E302        49 Male   10.4      12.259      0.641
## 20      F2402        47  Fem    8.9      10.674      0.577

$\blacksquare$

4. Tworzenie nowych zmiennych dzięki `mutate()`

Stwórz zmienną bill_ratio będącą stosunkiem długości dzioba do jego głębokości. Użyj do tego skróconego zestawu dane2, w którym pozostawisz tylko kolumny morfometryczne oraz ID matki (dam) oraz gniazdo (fosternest):

dane2 <- select(dane, dam, fosternest, tarsus, weight, starts_with("bill"))
head(mutate(dane2, bill_ratio = round(bill_length/bill_depth, 3)),
     20)

##        dam fosternest tarsus weight bill_length bill_depth bill_ratio
## 1  R187557      F2102  13.55    8.9      10.885      0.597     18.233
## 2  R187559      F1902  15.07   10.5      12.056      0.682     17.677
## 3  R187568       A602  14.99    9.9      12.000      0.673     17.831
## 4  R187518      A1302  14.69    9.7      11.733      0.639     18.362
## 5  R187528      A2602  14.46    9.7      11.565      0.631     18.328
## 6  R187945      C2302  13.93    9.3      11.165      0.626     17.835
## 7     Fem3      C1902  13.93    9.3      11.154      0.591     18.873
## 8  R187030      C1302  15.45   10.2      12.348      0.675     18.293
## 9  R187517       C602  14.31    9.5      11.434      0.629     18.178
## 10 R187523      B2202  14.46    9.6      11.564      0.609     18.989
## 11 R186902      B1402  14.99   10.0      11.977      0.636     18.832
## 12 R187400      B1002  14.08    9.7      11.251      0.685     16.425
## 13 R187932       B502  14.77    9.9      11.828      0.638     18.539
## 14 R187582      D1202  13.63    9.2      10.904      0.612     17.817
## 15 R187545      D1002  16.81   11.4      13.399      0.797     16.812
## 16 R187546       D902  14.24    9.8      11.409      0.679     16.803
## 17 R187590       D202  14.92   10.1      11.929      0.700     17.041
## 18 R187548       E902  14.46    9.6      11.512      0.669     17.208
## 19 R187594       E302  15.37   10.4      12.259      0.641     19.125
## 20 R187588      F2402  13.40    8.9      10.674      0.577     18.499

Zadanie Za pomocą funkcji paste() stwórz zmienną zawierającą kombinację ID matki (dam) i gniazda wychowywania (fosternest). Nazwij tą zmienną “crossfoster”.

Spodziewany wynik

head(mutate(dane2, crossfoster = paste(dam, fosternest, sep = "_")),
     20)

##        dam fosternest tarsus weight bill_length bill_depth   crossfoster
## 1  R187557      F2102  13.55    8.9      10.885      0.597 R187557_F2102
## 2  R187559      F1902  15.07   10.5      12.056      0.682 R187559_F1902
## 3  R187568       A602  14.99    9.9      12.000      0.673  R187568_A602
## 4  R187518      A1302  14.69    9.7      11.733      0.639 R187518_A1302
## 5  R187528      A2602  14.46    9.7      11.565      0.631 R187528_A2602
## 6  R187945      C2302  13.93    9.3      11.165      0.626 R187945_C2302
## 7     Fem3      C1902  13.93    9.3      11.154      0.591    Fem3_C1902
## 8  R187030      C1302  15.45   10.2      12.348      0.675 R187030_C1302
## 9  R187517       C602  14.31    9.5      11.434      0.629  R187517_C602
## 10 R187523      B2202  14.46    9.6      11.564      0.609 R187523_B2202
## 11 R186902      B1402  14.99   10.0      11.977      0.636 R186902_B1402
## 12 R187400      B1002  14.08    9.7      11.251      0.685 R187400_B1002
## 13 R187932       B502  14.77    9.9      11.828      0.638  R187932_B502
## 14 R187582      D1202  13.63    9.2      10.904      0.612 R187582_D1202
## 15 R187545      D1002  16.81   11.4      13.399      0.797 R187545_D1002
## 16 R187546       D902  14.24    9.8      11.409      0.679  R187546_D902
## 17 R187590       D202  14.92   10.1      11.929      0.700  R187590_D202
## 18 R187548       E902  14.46    9.6      11.512      0.669  R187548_E902
## 19 R187594       E302  15.37   10.4      12.259      0.641  R187594_E302
## 20 R187588      F2402  13.40    8.9      10.674      0.577 R187588_F2402

$\blacksquare$

5. Grupowanie obserwacji przez `group_by()`

Grupowanie pozwala stworzyć w zestawie danych “podzbiory” traktowane jako oddzielne pule przy liczeniu określonych podsumowań.

Porównaj taką operację:

head(mutate(dane2, mean_tars = mean(tarsus)),
     20)

##        dam fosternest tarsus weight bill_length bill_depth mean_tars
## 1  R187557      F2102  13.55    8.9      10.885      0.597   14.5005
## 2  R187559      F1902  15.07   10.5      12.056      0.682   14.5005
## 3  R187568       A602  14.99    9.9      12.000      0.673   14.5005
## 4  R187518      A1302  14.69    9.7      11.733      0.639   14.5005
## 5  R187528      A2602  14.46    9.7      11.565      0.631   14.5005
## 6  R187945      C2302  13.93    9.3      11.165      0.626   14.5005
## 7     Fem3      C1902  13.93    9.3      11.154      0.591   14.5005
## 8  R187030      C1302  15.45   10.2      12.348      0.675   14.5005
## 9  R187517       C602  14.31    9.5      11.434      0.629   14.5005
## 10 R187523      B2202  14.46    9.6      11.564      0.609   14.5005
## 11 R186902      B1402  14.99   10.0      11.977      0.636   14.5005
## 12 R187400      B1002  14.08    9.7      11.251      0.685   14.5005
## 13 R187932       B502  14.77    9.9      11.828      0.638   14.5005
## 14 R187582      D1202  13.63    9.2      10.904      0.612   14.5005
## 15 R187545      D1002  16.81   11.4      13.399      0.797   14.5005
## 16 R187546       D902  14.24    9.8      11.409      0.679   14.5005
## 17 R187590       D202  14.92   10.1      11.929      0.700   14.5005
## 18 R187548       E902  14.46    9.6      11.512      0.669   14.5005
## 19 R187594       E302  15.37   10.4      12.259      0.641   14.5005
## 20 R187588      F2402  13.40    8.9      10.674      0.577   14.5005

z taką:

dane3 <- group_by(dane2, dam)
head(mutate(dane3, mean_tars = mean(tarsus)),
     20)

## # A tibble: 20 × 7
## # Groups:   dam [20]
##    dam     fosternest tarsus weight bill_length bill_depth mean_tars
##    <fct>   <fct>       <dbl>  <dbl>       <dbl>      <dbl>     <dbl>
##  1 R187557 F2102        13.6    8.9        10.9      0.597      13.9
##  2 R187559 F1902        15.1   10.5        12.1      0.682      14.8
##  3 R187568 A602         15.0    9.9        12        0.673      14.3
##  4 R187518 A1302        14.7    9.7        11.7      0.639      14.5
##  5 R187528 A2602        14.5    9.7        11.6      0.631      14.4
##  6 R187945 C2302        13.9    9.3        11.2      0.626      14.7
##  7 Fem3    C1902        13.9    9.3        11.2      0.591      14.5
##  8 R187030 C1302        15.4   10.2        12.3      0.675      15.2
##  9 R187517 C602         14.3    9.5        11.4      0.629      14.8
## 10 R187523 B2202        14.5    9.6        11.6      0.609      14.2
## 11 R186902 B1402        15.0   10          12.0      0.636      14.2
## 12 R187400 B1002        14.1    9.7        11.3      0.685      14.7
## 13 R187932 B502         14.8    9.9        11.8      0.638      14.2
## 14 R187582 D1202        13.6    9.2        10.9      0.612      14.5
## 15 R187545 D1002        16.8   11.4        13.4      0.797      14.7
## 16 R187546 D902         14.2    9.8        11.4      0.679      14.6
## 17 R187590 D202         14.9   10.1        11.9      0.7        15.0
## 18 R187548 E902         14.5    9.6        11.5      0.669      14.2
## 19 R187594 E302         15.4   10.4        12.3      0.641      14.7
## 20 R187588 F2402        13.4    8.9        10.7      0.577      14.3

Aby lepiej zobaczyć co właściwie się tutaj stało - posortujmy wynikową tabelę wg ID matki:

dane3 <- group_by(dane2, dam)
dane3 <- mutate(dane3, mean_tars = mean(tarsus))
arrange(dane3, dam)

## # A tibble: 828 × 7
## # Groups:   dam [106]
##    dam   fosternest tarsus weight bill_length bill_depth mean_tars
##    <fct> <fct>       <dbl>  <dbl>       <dbl>      <dbl>     <dbl>
##  1 Fem2  H702         14.5   10          11.6      0.645      14.4
##  2 Fem2  H702         13.8    9.5        11.1      0.619      14.4
##  3 Fem2  A18B02       14.4    9.5        11.5      0.634      14.4
##  4 Fem2  A18B02       14.7    9.5        11.8      0.658      14.4
##  5 Fem2  A18B02       14.8    9.8        11.8      0.607      14.4
##  6 Fem20 H1102        13.9    9.4        11.1      0.567      14.5
##  7 Fem20 H1102        14.5    9.9        11.6      0.665      14.5
##  8 Fem20 H1102        14.8    9.9        11.9      0.643      14.5
##  9 Fem20 H1102        14.1    9.4        11.2      0.591      14.5
## 10 Fem20 H1102        14.6    9.8        11.7      0.643      14.5
## # … with 818 more rows

Czy zauważyłeś zmianę tabeli danych na nowy typ obiektu zwany tibble?

Wykonanie powyższej serii operacji bardzo ułatwia użycie tzw. operatora pipe, pozwalającego łączyć dłuższe ciągi operacji w jeden spójny “potok”:

group_by(dane2, dam) %>%
  mutate(mean_tars = mean(tarsus)) %>%
  arrange(dam)

## # A tibble: 828 × 7
## # Groups:   dam [106]
##    dam   fosternest tarsus weight bill_length bill_depth mean_tars
##    <fct> <fct>       <dbl>  <dbl>       <dbl>      <dbl>     <dbl>
##  1 Fem2  H702         14.5   10          11.6      0.645      14.4
##  2 Fem2  H702         13.8    9.5        11.1      0.619      14.4
##  3 Fem2  A18B02       14.4    9.5        11.5      0.634      14.4
##  4 Fem2  A18B02       14.7    9.5        11.8      0.658      14.4
##  5 Fem2  A18B02       14.8    9.8        11.8      0.607      14.4
##  6 Fem20 H1102        13.9    9.4        11.1      0.567      14.5
##  7 Fem20 H1102        14.5    9.9        11.6      0.665      14.5
##  8 Fem20 H1102        14.8    9.9        11.9      0.643      14.5
##  9 Fem20 H1102        14.1    9.4        11.2      0.591      14.5
## 10 Fem20 H1102        14.6    9.8        11.7      0.643      14.5
## # … with 818 more rows

Zadanie Korzystając z operatora %>% zmodyfikuj powyższy ciąg tak, by na końcu powstała tabela o następujących właściwościach:

grupowanie ze względu na fosternest;
obliczenie średniej wagi piskląt w obrębie każdego gniazda wychowywania;
usunięcie z tabeli zmiennych niepotrzebnych (cechy dzioba)
obliczenie stosunku wagi średniej w danym gnieździe do średniej wagi w całej populacji. (Uwaga! Ten etap wymaga usunięcia grupowania za pomocą funkcji ungroup())

Aby wyświetlić więcej niż 10 wierszy tabeli dodaj do potoku funkcji print() z opcją n = 30 (jeśli chcesz wyświetlić 30 wierszy). Dlaczego w trakcie obliczeń musimy zastosować ungroup()? Co się stanie, jeśli tego nie zrobimy?

Spodziewany wynik

group_by(dane2, fosternest) %>%
  mutate(mean_weight = mean(weight)) %>%
  select(-starts_with("bill")) %>%
  arrange(fosternest) %>%
  ungroup() %>%
  mutate(mean_ratio = mean_weight/mean(weight)) %>%
  print(n = 30)

## # A tibble: 828 × 6
##    dam     fosternest tarsus weight mean_weight mean_ratio
##    <fct>   <fct>       <dbl>  <dbl>       <dbl>      <dbl>
##  1 R187537 A1002        14.7   10          9.64      0.985
##  2 R187537 A1002        14.8    9.8        9.64      0.985
##  3 R187537 A1002        14.9    9.7        9.64      0.985
##  4 R187537 A1002        14.5    9.6        9.64      0.985
##  5 R187563 A1002        13.2    9.3        9.64      0.985
##  6 R187563 A1002        13.6    9.2        9.64      0.985
##  7 R187563 A1002        14.0    9.2        9.64      0.985
##  8 R187563 A1002        14.7   10.1        9.64      0.985
##  9 R187563 A1002        14.4    9.9        9.64      0.985
## 10 R187539 A102         14.2    9.8        9.67      0.987
## 11 R187539 A102         14.1    9.5        9.67      0.987
## 12 R187539 A102         15.1   10.1        9.67      0.987
## 13 R187539 A102         13.1    8.8        9.67      0.987
## 14 R187566 A102         14.6    9.4        9.67      0.987
## 15 R187566 A102         14.7    9.9        9.67      0.987
## 16 R187566 A102         14.4    9.6        9.67      0.987
## 17 R187566 A102         14.5    9.8        9.67      0.987
## 18 R187566 A102         15.1   10.1        9.67      0.987
## 19 R186910 A1202        14.5   10.2        9.7       0.991
## 20 R186910 A1202        13.5    9.1        9.7       0.991
## 21 R186910 A1202        14.8   10.4        9.7       0.991
## 22 R186910 A1202        14.1    9.2        9.7       0.991
## 23 R186910 A1202        15.1   10.3        9.7       0.991
## 24 R186910 A1202        14.1    9.5        9.7       0.991
## 25 R186910 A1202        13.9    9.2        9.7       0.991
## 26 R187518 A1302        14.7    9.7        9.84      1.00 
## 27 R187518 A1302        14.8   10.1        9.84      1.00 
## 28 R187518 A1302        14.1    9.4        9.84      1.00 
## 29 R187518 A1302        14.6    9.9        9.84      1.00 
## 30 R187518 A1302        14.1    9.7        9.84      1.00 
## # … with 798 more rows

$\blacksquare$

6. Podsumowanie z `summarise()`

Przefiltruj dane tak, by zostały w nich tylko ptaki ze znaną płcią (sex == "Male" lub sex == "Fem") i podsumuj je wyliczając dla obserwacji średnią długość skoku, liczbę obserwacji oraz błąd standardowy skoku (liczony jako SD/sqrt(N)):

dane %>%
  filter(sex != "UNK") %>%
  summarise(mean_t = mean(tarsus),
            N = n(),
            se_t = sd(tarsus)/sqrt(n()))

##     mean_t   N       se_t
## 1 14.50553 781 0.01797891

Zadanie Powtórz powyższe obliczenia, ale dla ptaków w każdej możliwej kombinacji sex oraz habitat (tzn. samce w lesie, samce w parku, etc.)

Spodziewany wynik

dane %>%
  filter(sex != "UNK") %>%
  group_by(sex, habitat) %>%
  summarise(mean_t = mean(tarsus),
            N = n(),
            se_t = sd(tarsus)/sqrt(n()))

## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.

## # A tibble: 4 × 5
## # Groups:   sex [2]
##   sex   habitat mean_t     N   se_t
##   <fct> <fct>    <dbl> <int>  <dbl>
## 1 Fem   forest    14.3   188 0.0334
## 2 Fem   park      14.3   185 0.0343
## 3 Male  forest    14.7   205 0.0325
## 4 Male  park      14.7   203 0.0325

$\blacksquare$

7. Szerokie i długie zestawy danych

Oryginalny zestaw dane jest typu “długiego”: każda obserwacja to jeden osobnik, a zmienne kategoryczne zajmują po jednej kolumnie zawierającej przypisania każdego osobnika do określonej kategorii. Często jednak dane, które dostajemy są typu “szerokiego” (np. grupy eksperymentalne czy płcie znajdują się w osobnych kolumnach). Między tymi dwoma typami danych można sprawnie się przemieszczać za pomocą funkcji pivot_wider() i pivot_longer().

By dowiedzieć się o nich nieco więcej - zajrzyj tutaj:

vignette("pivot")

Stwórzmy przykładowy zestaw danych zawierający średnie ciężary piskląt i ich wariancje według płci, pogrupowane według fosternest:

dane_l <- dane %>%
  filter(sex != "UNK") %>%
  group_by(fosternest, sex) %>%
  summarise(mean_w = mean(weight), var_w = var(weight))

## `summarise()` has grouped output by 'fosternest'. You can override using the
## `.groups` argument.

Taki zestaw danych w formie długiej możemy przerobić na szeroki tworząc osobne kolumny dla samców i samic w danym gnieździe.

dane_l %>%
  pivot_wider(names_from = sex, values_from = mean_w)

## # A tibble: 201 × 4
## # Groups:   fosternest [104]
##    fosternest   var_w   Fem  Male
##    <fct>        <dbl> <dbl> <dbl>
##  1 A1002      0.139    9.57 NA   
##  2 A1002      0.07    NA     9.8 
##  3 A102       0.249    9.48 NA   
##  4 A102       0.0770  NA     9.82
##  5 A1202      0.0433   9.27 NA   
##  6 A1202      0.00500 NA    10.4 
##  7 A1302      0        9.7  NA   
##  8 A1302      0.0817  NA     9.88
##  9 A1602      0.0960   9.5  NA   
## 10 A1602      0.180   NA     9.7 
## # … with 191 more rows

Powyższa komenda nie zadziałała tak, jak chcemy, bo pozostawiliśmy w danych kolumnę z wariancją, zawierającą unikatowe wpisy “dublujące” wiersze z każdym fosternest. Musimy albo się tej kolumny pozbyć, albo też uwzględnić ją w pivot_wider:

dane_l %>%
  select(-var_w) %>%
  pivot_wider(names_from = sex, values_from = mean_w)

## # A tibble: 104 × 3
## # Groups:   fosternest [104]
##    fosternest   Fem  Male
##    <fct>      <dbl> <dbl>
##  1 A1002       9.57  9.8 
##  2 A102        9.48  9.82
##  3 A1202       9.27 10.4 
##  4 A1302       9.7   9.88
##  5 A1602       9.5   9.7 
##  6 A1802       9.6   9.95
##  7 A18B02     10     9.72
##  8 A2202       9.48 10   
##  9 A22B02      9.33  9.48
## 10 A2302       9.57  9.7 
## # … with 94 more rows

dane_l %>%
  pivot_wider(names_from = sex, values_from = c(mean_w, var_w))

## # A tibble: 104 × 5
## # Groups:   fosternest [104]
##    fosternest mean_w_Fem mean_w_Male var_w_Fem var_w_Male
##    <fct>           <dbl>       <dbl>     <dbl>      <dbl>
##  1 A1002            9.57        9.8     0.139     0.07   
##  2 A102             9.48        9.82    0.249     0.0770 
##  3 A1202            9.27       10.4     0.0433    0.00500
##  4 A1302            9.7         9.88    0         0.0817 
##  5 A1602            9.5         9.7     0.0960    0.180  
##  6 A1802            9.6         9.95    0.260     0.605  
##  7 A18B02          10           9.72    0.0200    0.0870 
##  8 A2202            9.48       10       0.192     0.167  
##  9 A22B02           9.33        9.48    0.0233    0.0692 
## 10 A2302            9.57        9.7     0.171     0.07   
## # … with 94 more rows

Zapiszmy dane w nowym obiekcie dane_w.

dane_w <- dane_l %>%
  select(-var_w) %>%
  pivot_wider(names_from = sex, values_from = mean_w)

Zadanie Spróbuj odwrócić działanie zamieniając “szeroki” styl danych na “długi” - korzystając z funkcji pivot_longer(). Nazwij zmienną “płeć” animal_sex, a zmienną z ciężarem ciała - ’body_weight`.

Spodziewany wynik

dane_w %>%
  pivot_longer(cols = c(Fem, Male),
               names_to = "animal_sex",
               values_to = "body_weight")

## # A tibble: 208 × 3
## # Groups:   fosternest [104]
##    fosternest animal_sex body_weight
##    <fct>      <chr>            <dbl>
##  1 A1002      Fem               9.57
##  2 A1002      Male              9.8 
##  3 A102       Fem               9.48
##  4 A102       Male              9.82
##  5 A1202      Fem               9.27
##  6 A1202      Male             10.4 
##  7 A1302      Fem               9.7 
##  8 A1302      Male              9.88
##  9 A1602      Fem               9.5 
## 10 A1602      Male              9.7 
## # … with 198 more rows

B. Prosta statystyka

1. Czy samce są większe od samic?

W R proste testy statystyczne wykonujemy za pomocą serii funkcji z rodziny __.test - np. takich:

Funkcja	Test	Zastosowanie
`t.test()`	test t-Studenta	porównywanie średnich
`var.test()`	test F	porównywanie wariancji (zmienności)
`cor.test()`	test korelacji	testowanie siły związku między zmiennymi
`chisq.test()`	test $\chi^2$	testowanie zmiennych kategorycznych i tabel wielopolowych
`kruskal.test()`	test Kruskall-Wallisa	nieparametryczne porównywanie grup
`wilcox.test`	test Wilcoxona	nieparametryczny test par wiązanych

Sprawdźmy za pomocą testu korelacji czy waga piskląt związana jest istotnie z ich długością skoku. Jaka jest: wartość p testu, obliczona korelacja i jej 95% przedział ufności? Czy odrzucasz hipotezę zerową? Jaka ona jest?

test1 <- cor.test(dane$tarsus, dane$weight)
test1

## 
##  Pearson's product-moment correlation
## 
## data:  dane$tarsus and dane$weight
## t = 48.044, df = 826, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8391000 0.8751343
## sample estimates:
##       cor 
## 0.8581706

test1$p.value

## [1] 2.096071e-241

Prosty wykres rozrzutu dobrze obrazuje otrzymany wynik:

plot(dane$tarsus, dane$weight,
     xlab = "Weight", ylab = "Tarsus")

Zadanie Wiedząc, że opcja wykresu col może m.in. przyjmować wartości będące zmiennymi kategorycznymi - jak spróbujesz dodać do wykresy informację o płci, zmieniającą kolor obserwacji w zależności od płci? Wskazówka Zamiast dane$zmienna - użyj jitter(dane$zmienna, 3).

Sugerowany wynik

plot(jitter(dane$tarsus, 3), jitter(dane$weight, 3),
     xlab = "Weight", ylab = "Tarsus",
     col = dane$sex)

Zadanie Z wykresu jasno wynika, że samce powinny być większe od samic pod względem zarówno masy, jak i skoku. Przetestujmy tą różnicę za pomocą testu t-Studenta. Wykonaj go ekstrahując z danych skok samic, a potem samców (za pomocą funkcji subset(), lub w inny znany Ci sposób). Dla przypomnienia: subset(dane, habitat == "forest")$weight pozwoliłby wyciągnąć z danych wagi ptaków żyjących w lesie. Wykonaj go w podobny sposób, jak cor.test(). Jaka jest Twoja decyzja: odrzucasz czy przyjmujesz hipotezę zerową o braku różnic między płciami?

Sugerowany wynik

test2 <- t.test(subset(dane, sex == "Male")$tarsus,
                subset(dane, sex == "Fem")$tarsus)
test2

## 
##  Welch Two Sample t-test
## 
## data:  subset(dane, sex == "Male")$tarsus and subset(dane, sex == "Fem")$tarsus
## t = 11.759, df = 773.68, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.3250593 0.4553372
## sample estimates:
## mean of x mean of y 
##  14.69189  14.30169

W przypadku danych w formacie “długim” posługiwanie się funkcją t.test() jest znacznie łatwiejsze w konwencji tzw. formuły. Formuła to rodzaj wyrażenia matematycznego, które pozwala zapisać zależności istniejące między zmienną/zmiennymi zależną/zależnymi i niezależną/niezależnymi. W przypadku prostej zależności, gdzie mamy jedną zmienną zależną y (np. skok) i jedną zmienną niezależną x (np. płeć) - formułę możemy zapisać jako y ~ x.

Zadanie Powtórz powyższy test t-Studenta używają notacji formuły zamiast dwóch osobnych wektorów zawierających dane samic i samców. Aby konwencja ta zadziałała będziesz potrzebować wyrażenia opisującego zależność skoku od płci (tarsus ~ sex) oraz dodatkowej opcji data w t.test() wskazującej na to z jakiego obiektu pochodzą dane (zajrzyj do ?t.test() w razie wątpliwości). Uwaga! W dalszym ciągu musisz usunąć z danych obserwacje sex == "UNK" - w teście t-Studenta możemy porównywać tylko 2 grupy.

Sugerowany wynik

test3 <- t.test(tarsus ~ sex, data = subset(dane, sex != "UNK"))
test3

## 
##  Welch Two Sample t-test
## 
## data:  tarsus by sex
## t = -11.759, df = 773.68, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Fem and group Male is not equal to 0
## 95 percent confidence interval:
##  -0.4553372 -0.3250593
## sample estimates:
##  mean in group Fem mean in group Male 
##           14.30169           14.69189

Zadanie Przeprowadzany test t-Studenta zakłada, że wariancje porównywanych grupo są różne. Przeprowadzając test funkcją var.test() zweryfikuj to założenie, testując hipotezę zerową, że wariancje są równe. Składnia testu jest bardzo podobna do tej użytej powyżej z wykorzystaniem formuły.

Sugerowany wynik

test4 <- var.test(tarsus ~ sex, data = subset(dane, sex != "UNK"))
test4

## 
##  F test to compare two variances
## 
## data:  tarsus by sex
## F = 0.98626, num df = 372, denom df = 407, p-value = 0.8928
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.808562 1.204097
## sample estimates:
## ratio of variances 
##          0.9862598

$\blacksquare$

2. Czy samic lub samców jest więcej w lesie/w parku?

Aby odpowiedzieć na to pytanie - potrzebujemy podsumowania zliczającego samce i samice w dwóch rodzajach siedliska. Najłatwiej uzyskać takie podsumowanie stosując funkcję table():

dane2 <- filter(dane, sex != "UNK") %>% gdata::drop.levels()
sex_hab <- table(dane2$sex, dane2$habitat)
sex_hab

##       
##        forest park
##   Fem     188  185
##   Male    205  203

Zadanie Wiedząc, że funkcja chisq.test potrafi jako jeden z argumentów przyjmować całą tabelę wielopolową - przeprowadź test chi-kwadrat hipotezy badającej zależność występowania ptaków w 2 siedliskach w zależności od płci. Czy zależność ta jest istotna? Co znajduje się w zbiorze stdres powstałego obiektu?

Spodziewany wynik

test5 <- chisq.test(sex_hab)
test5

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  sex_hab
## X-squared = 0, df = 1, p-value = 1

test5$stdres

##       
##             forest        park
##   Fem   0.04384568 -0.04384568
##   Male -0.04384568  0.04384568

$\blacksquare$

3. Czy długość skoku zależy od daty wyklucia? Regresja prosta

Sprawdźmy to najpierw wizualnie:

plot(tarsus ~ hatchdate, data = dane)

Nie wygląda to na zbyt silną zależność - możemy jednak sprawdzić to formalnie za pomocą modelu liniowego w najprostszej jego wersji. Model liniowy, jaki będziemy testować wygląda następująco:

\[y_i = a + bx_i + e_i\] czyli

\[tarsus_i = a + b*hatchdate_i + e_i\] Do jego obliczenia służy funkcja lm() (linear model) - co zaskakujące przyjmuje ona argumenty w sposób bardzo podobny do np. t.test() (korzystając z konwencji formuły).

Zadanie Stosując analogię wywołania funkcji t.test(), użyj lm() by przetestować zależność długości skoku sikory modrej od daty klucia (tarsus ~ hatchdate). W pierwszym kroku stwórz obiekt model zawierający wynik działania funkcji lm(). W drugim kroku użyj summary() by podejrzeć wynik działania obliczeń.

Sugerowany wynik

model <- lm(tarsus ~ hatchdate, data = dane)
summary(model)

## 
## Call:
## lm(formula = tarsus ~ hatchdate, data = dane)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8746 -0.3471  0.0468  0.3011  2.2925 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 15.08869    0.33443  45.117   <2e-16 ***
## hatchdate   -0.01215    0.00690  -1.761   0.0786 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4996 on 826 degrees of freedom
## Multiple R-squared:  0.003741,   Adjusted R-squared:  0.002535 
## F-statistic: 3.102 on 1 and 826 DF,  p-value: 0.07858

Przyjrzyj się wykresom, które powstają dzięki wywołaniu funkcji

plot(model)

Pozwalają one ocenić, czy nasz model statystyczny spełnia założenia regresji prostej. Który z wykresów najlepiej oddaje test wizualny założenia zakładającego normalność wartości resztkowych modelu?

$\blacksquare$

C. Expressem przez bootstrapping

Do tej pory do naszych analiz wykorzystywaliśmy metody parametryczne - możemy jednak spróbować wykorzystać inną, prostszą metodę, która pozwoli nam uwolnić się od ograniczeń rozkładowych naszych danych. Spójrzmy na przykład na rozkład zmiennej hatchdate - nie jest on najpiękniejszy i na pierwszy rzut oka nie przypomina on zgrabnego rozkładu normalnego:

hist(dane$hatchdate, 15)

Spróbujmy skonstruować test oparty na bootstrappingu, którym sprawdzimy czy ptaki z lasu (forest) różnią się pod kątem daty klucia od ptaków z parków (park). Szybki rzut oka na wykres pudełkowy nie wskazuje na zbyt wielkie różnice.

plot(hatchdate ~ habitat, data = dane)

Do skonstruowania bootstrappingu wykorzystamy pętlę for, która pozwala powtarzać jakąś operację wiele razy. Jeśli chcielibyśmy np. wydrukować na ekranie (print()) jakiś tekst pięć razy - moglibyśmy osiągnąć to w następujący sposób:

for(i in 1:5) {
  print("Tekst najbardziej oryginalny!")
}

## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"

Zauważ, że pętla for wykorzystuje jakąś zmienną (tutaj i), zmieniając jej wartość w zadanym zakresie (tutaj 1:5) i za każdym razem (dla każ∂ego i) wykonując jakieś działania ujęte wewnątrz {}.

Potrzebna nam będzie również metoda losowego wybierania pewnych wierszy z naszych oryginalnych danych - ze zwracaniem. Idealnie do tego nadaje się funkcja sample(). Zobaczmy, jak losuje ona zadaną liczbę wartości (10) z podanej listy (liczby 1 do 20), ze zwracaniem (replace = T):

sample(1:20, size = 10, replace = T)

##  [1] 13  1  2 15 10 16 19 13  1 15

Uwaga! Funkcja ta używa generatora liczb pseudolosowych - więc u każdego wynik tej operacji będzie inny!

Wreszcie - musimy mieć naszą statystykę testową. Niech będzie nią różnica między średnią datą klucia w sikor leśnych i parkowych $d_0$ - jeśli wynosi ona zero, ptaki mają identyczne daty klucia:

d0 <- mean(filter(dane, habitat == "forest")$hatchdate) -
  mean(filter(dane, habitat == "park")$hatchdate)
d0

## [1] -0.1695764

Ujemna różnica wskazuje na nieco wyższą wartość (późniejsze klucie) w parku - pytanie, czy istotnie późniejsze?

Oto co powinna robić nasza procedura bootstrappingu:

Krok 1: wylosuj z zestawu danych k = 828 obserwacji (wierszy) ze zwracaniem (tyle wierszy składa się na oryginalny zestaw danych).
Krok 2: dla każdego z tych wylosowanych zestawów policz jego własną statystykę d - różnicę w średniej dacie klucia między lasem i parkiem.
Krok 3: Powtórz całą procedurę N = 1000 razy (jest to arbitralna wielkość, im więcej - tym lepiej). Za każdym razem zapisz wynik do zbiorczego obiekty out.
Krok 4: Skonstruuj rozkład statystyki d i zobacz gdzie wypada w niej jej oryginalna (oparta na wejściowych danych) wartość d0.

N <- 1000 # liczba zakładanych powtórzeń
k <- nrow(dane) # wielkość oryginalnego zbioru danych

# oryginalna statystyka testowa
d0 <- mean(filter(dane, habitat == "forest")$hatchdate) -
  mean(filter(dane, habitat == "park")$hatchdate)

# obiekt wynikowy przygotowany na N=1000 końcowych wartości
out <- numeric(N)

for(i in 1:N) {
  
  # losujemy 828 wierszy, ze zwracaniem
  # a oryginalnej puli unikatowych 828 wierszy danych
  rows <- sample(1:k, size = k, replace = T)
  
  # tworzymy tymczasowe dane wewnątrz pętli
  # zawierające tylko wylosowane wiersze
  temp_dane <- dane[rows, ]
  
  # oblicza na podstawie tak przygotowanych danych
  # tymczasową statystykę testową dla powtórzenia
  # nr i pętli
  d_temp <- mean(filter(temp_dane, habitat == "forest")$hatchdate) -
    mean(filter(temp_dane, habitat == "park")$hatchdate)
  
  # zapisujemy wyliczoną statystykę
  # do zmiennej out w miejscu i
  out[i] <- d_temp
  
}

Wyświetlmy rozkład symulowanej statystyki d oraz położenie d0 na rozkładzie:

hist(out, 50)
abline(v = d0, lwd = 2)
abline(v = 0, col = "red", lwd = 2, lty = 2)

Jak podejdziesz to odpowiedzi na pytanie - czy ta symulacja wskazuje na istotną statystycznie różnicę między symulowaną statystyką a wartością zero?

Class 2

Zanim zaczniemy

0. Przydatne zasoby pomocowe w pracy z R - część druga

A. tidyverse oraz dplyr - “porządne dane”

1. Filtrowanie danych: filter()

2. Sortowanie za pomocą arrange()

3. Selekcja zmiennych dzięki select()

4. Tworzenie nowych zmiennych dzięki mutate()

5. Grupowanie obserwacji przez group_by()

6. Podsumowanie z summarise()

7. Szerokie i długie zestawy danych

B. Prosta statystyka

1. Czy samce są większe od samic?

2. Czy samic lub samców jest więcej w lesie/w parku?

3. Czy długość skoku zależy od daty wyklucia? Regresja prosta

C. Expressem przez bootstrapping

0. Przydatne zasoby pomocowe w pracy z `R` - część druga

A. `tidyverse` oraz `dplyr` - “porządne dane”

1. Filtrowanie danych: `filter()`

2. Sortowanie za pomocą `arrange()`

3. Selekcja zmiennych dzięki `select()`

4. Tworzenie nowych zmiennych dzięki `mutate()`

5. Grupowanie obserwacji przez `group_by()`

6. Podsumowanie z `summarise()`